Load all required libraries
If one or more packages is not already installed for you, you can
install it using install.packages("PackageName")
library(readr)
library(dplyr)
library(stringr)
library(survival)
library(ggplot2)
library(survminer)
library(wesanderson)
library(ggpubr)
library(coxme)
library(lme4)
library(corrplot)
library(RColorBrewer)
library(ggrepel)
library(reshape2)
library(ggcorrplot)
library(nlme)
Set your working directory and import dataset
Set your working directory to the folder where your dataset is
present
The dataset can be downloaded from Main
Data and Bacterial
Load Data
setwd("~/Documents/GitHub/Pathogen-Potential")
Data <- read_csv("PathogenPotential_AG_04062023.csv")
LoadData <- read_csv("PathogenPotential_BacterialLoad_AG.csv")
Data processing and transformation to survival format
Entering Dpt host genotype information
Data <- Data %>%
mutate(Genotype = recode(Genotype,
"1G" = "DptNull",
"1N" = "DptR",
"8N" = "DptS"))
Converting survival data to 0/1 format
#Create an empty dataframe Data_1 to store the survival information
Data_1 <- data.frame(Vial_ID = character(),
Genotype = character(),
Treatment = character(),
OD_group = character(),
Date = character(),
Day_Dead = numeric(),
censor = numeric())
#Convert data to 0/1 format
for(i in 1:length(Data$Vial_ID)){
temp=data.frame(Data[i,])
for(j in seq(10,22,2)){
k = temp[1,j-2]-temp[1,j]
if(k>0){
for(kk in 1:k){
Data_1[length(Data_1$Vial_ID)+1,]=c(temp[1,1],
as.character(temp[1,2]),
as.character(temp[1,3]),
temp[1,5],
as.character(temp[1,6]),
j/2 - 4,
1)
}
}
}
for(j in 8){
k = temp[1,j]
if(k>0){
for(kk in 1:k){
Data_1[length(Data_1$Vial_ID)+1,]=c(temp[1,1],
as.character(temp[1,2]),
as.character(temp[1,3]),
temp[1,5],
as.character(temp[1,6]),
j/2 - 4,
0)
}
}
}
for(j in 22){
k = temp[1,j]
if(k>0){
for(kk in 1:k){
Data_1[length(Data_1$Vial_ID)+1,]=c(temp[1,1],
as.character(temp[1,2]),
as.character(temp[1,3]),
temp[1,5],
as.character(temp[1,6]),
j/2 - 4,
0)
}
}
}
}
#Converting Data_1 columns from character to numeric/factors
Data_1$Day_Dead = as.numeric(Data_1$Day_Dead)
Data_1$censor = as.numeric(Data_1$censor)
Data_1$Genotype <- as.factor(Data_1$Genotype)
Data_1$Treatment <- as.factor(Data_1$Treatment)
Data_1$Date <- as.factor(Data_1$Date)
Data_1$OD_group <- as.factor(Data_1$OD_group)
#Data_1 is the dataframe containing only survival data
Calculating fraction symptomatic
Dead individuals on any day are assumed to be symptomatic
Data$Fs_Day1 <- (Data$Day0_Alive - Data$Day1_Climbing) / Data$Day0_Alive
Data$Fs_Day2 <- (Data$Day0_Alive - Data$Day2_Climbing) / Data$Day0_Alive
Data$Fs_Day3 <- (Data$Day0_Alive - Data$Day3_Climbing) / Data$Day0_Alive
Data$Fs_Day4 <- (Data$Day0_Alive - Data$Day4_Climbing) / Data$Day0_Alive
Data$Fs_Day5 <- (Data$Day0_Alive - Data$Day5_Climbing) / Data$Day0_Alive
Data$Fs_Day6 <- (Data$Day0_Alive - Data$Day6_Climbing) / Data$Day0_Alive
Data$Fs_Day7 <- (Data$Day0_Alive - Data$Day7_Climbing) / Data$Day0_Alive
#Creating a dataframe Fs for fraction symptomatic by transforming data
Fs <- melt(Data[,c(1:3,5,6,24:30)], id=c("Vial_ID",
"Genotype",
"Treatment",
"OD_group",
"Date"))
colnames(Fs)[7] <- "Fs"
Fs[c('F', 'Day')] <- str_split_fixed(Fs$variable, '_Day', 2)
Fs <- Fs[,c(1:5,7,9)]
Calculating Pathogen Potential
PP = (Fs/I) * (10^M)
Data$PP_Day1 <- (Data$Fs_Day1/exp(Data$OD_group)) *
(10 ^ ((Data$Day0_Alive - Data$Day1_Alive)/Data$Day0_Alive))
Data$PP_Day2 <- (Data$Fs_Day2/exp(Data$OD_group)) *
(10 ^ ((Data$Day0_Alive - Data$Day2_Alive)/Data$Day0_Alive))
Data$PP_Day3 <- (Data$Fs_Day3/exp(Data$OD_group)) *
(10 ^ ((Data$Day0_Alive - Data$Day3_Alive)/Data$Day0_Alive))
Data$PP_Day4 <- (Data$Fs_Day4/exp(Data$OD_group)) *
(10 ^ ((Data$Day0_Alive - Data$Day4_Alive)/Data$Day0_Alive))
Data$PP_Day5 <- (Data$Fs_Day5/exp(Data$OD_group)) *
(10 ^ ((Data$Day0_Alive - Data$Day5_Alive)/Data$Day0_Alive))
Data$PP_Day6 <- (Data$Fs_Day6/exp(Data$OD_group)) *
(10 ^ ((Data$Day0_Alive - Data$Day6_Alive)/Data$Day0_Alive))
Data$PP_Day7 <- (Data$Fs_Day7/exp(Data$OD_group)) *
(10 ^ ((Data$Day0_Alive - Data$Day7_Alive)/Data$Day0_Alive))
#Creating a dataframe PP for pathogen potential by transforming data
PP <- melt(Data[,c(1:3,5,6,31:37)], id=c("Vial_ID",
"Genotype",
"Treatment",
"OD_group",
"Date"))
colnames(PP)[7] <- "PP"
PP[c('p', 'Day')] <- str_split_fixed(PP$variable, '_Day', 2)
PP <- PP[,c(1:5,7,9)]
Calculating proportion survival
#Proportion survival on Day X
Data$Prop_Day0_Alive <- Data$Day0_Alive/Data$Day0_Alive
Data$Prop_Day1_Alive <- Data$Day1_Alive/Data$Day0_Alive
Data$Prop_Day2_Alive <- Data$Day2_Alive/Data$Day0_Alive
Data$Prop_Day3_Alive <- Data$Day3_Alive/Data$Day0_Alive
Data$Prop_Day4_Alive <- Data$Day4_Alive/Data$Day0_Alive
Data$Prop_Day5_Alive <- Data$Day5_Alive/Data$Day0_Alive
Data$Prop_Day6_Alive <- Data$Day6_Alive/Data$Day0_Alive
Data$Prop_Day7_Alive <- Data$Day7_Alive/Data$Day0_Alive
Calculating cumulative risk scores https://stat.ethz.ch/R-manual/R-devel/library/survival/html/predict.coxph.html
Data_Risk <- subset(Data_1, !Day_Dead==0)
cox_model <- coxph(Surv(Day_Dead, censor) ~ as.factor(Vial_ID) +
Genotype +
Treatment +
OD_group,
data = Data_Risk)
# Compute the risk scores
risk_scores <- predict(cox_model, type = 'risk')
#Add it to a dataframe
Data_Risk$Risk_Score <- risk_scores
#group by vial id and day to assign cumulative scores
Data_1_2 <- Data_Risk %>%
group_by(Vial_ID,Genotype,Treatment,OD_group) %>%
summarise(Vial_ID=Vial_ID,
Genotype=Genotype,
Treatment=Treatment,
OD_group=OD_group,
Risk_Score=sum(Risk_Score)) %>%
distinct()
#Add risk scores to the main dataframe
Data <- merge(Data, Data_1_2)
Calculating PP_T & Fs_T & Fs/T
T=median survival day
#PP_T Fs_T Fs/T
#Estimating mean survival day
Surv <- melt(Data[,c(1:4,6,39:45)], id=c("Vial_ID",
"Genotype",
"Treatment",
"OD_group",
"Date"))
colnames(Surv)[7] <- "Prop_alive"
Surv[c('p', 'Day')] <- str_split_fixed(Surv$variable, '_Alive', 2)
Surv[c('p', 'Day')] <- str_split_fixed(Surv$p, 'Prop_Day', 2)
Surv <- Surv[,c(1:5,7,9)]
Surv$Day_T <- 7
for(i in 1:length(Surv$Vial_ID)){
temp=data.frame(Surv[i,])
for(j in 6) {
k=temp[1,j]
if(k<=0.5) {
Surv$Day_T[i] = Surv$Day[i]
}
}
}
Surv_1 <- Surv %>%
group_by(Vial_ID,Genotype,Treatment,OD_group) %>%
summarise(Vial_ID=Vial_ID,
Genotype=Genotype,
Treatment=Treatment,
OD_group=OD_group,
Day_T=min(Day_T)) %>%
distinct()
Data <- merge(Data, Surv_1)
Data$Fs_T <- NA
Data$PP_T <- NA
#Estimating pathogen potential and fraction symptomatic on mean survival day T
for(i in 1:length(Data$Vial_ID)) {
temp=data.frame(Data[i,])
j=temp[1,47]
if(j==1) {
Data$Fs_T[i] = Data$Fs_Day1[i]
Data$PP_T[i] = Data$PP_Day1[i]
}
if(j==2) {
Data$Fs_T[i] = Data$Fs_Day2[i]
Data$PP_T[i] = Data$PP_Day2[i]
}
if(j==3) {
Data$Fs_T[i] = Data$Fs_Day3[i]
Data$PP_T[i] = Data$PP_Day3[i]
}
if(j==4) {
Data$Fs_T[i] = Data$Fs_Day4[i]
Data$PP_T[i] = Data$PP_Day4[i]
}
if(j==5) {
Data$Fs_T[i] = Data$Fs_Day5[i]
Data$PP_T[i] = Data$PP_Day5[i]
}
if(j==6) {
Data$Fs_T[i] = Data$Fs_Day6[i]
Data$PP_T[i] = Data$PP_Day6[i]
}
if(j==7) {
Data$Fs_T[i] = Data$Fs_Day7[i]
Data$PP_T[i] = Data$PP_Day7[i]
}
}
Data$`Fs_T/T` <- Data$Fs_T/as.numeric(Data$Day_T)
Main Figures
Figure 1B
Fig1b <- ggplot(Data, aes(x = as.factor(OD_group),
y = Prop_Day4_Alive,
pch = as.factor(OD_group),
col = Genotype,
#fill = as.numeric(OD_group)
)) +
stat_summary(geom = "errorbar",
position = position_dodge(0.7),
width = 0.3,
#aes(alpha = as.numeric(OD_group))
) +
stat_summary(geom = "point",
position = position_dodge(0.7),
size = 6,
#aes(alpha = as.numeric(OD_group))
) +
facet_grid(~Treatment, scales = "free_x") +
coord_cartesian(ylim=c(0,1)) +
scale_x_discrete(
labels = function(x) {
if ("LB" %in% unique(Data$Treatment)) {
# For LB, show only '0'
ifelse(x == "0", "0", x)
} else {
# For other treatments, show the specific values
ifelse(x %in% c("0.01", "0.1", "0.5", "1", "4", "8"), x, NA)
}
}
) +
scale_shape_manual(values = c(4,1,2,0,5,7,8)) +
theme_classic() +
labs(x = "OD",
y = "Proportion alive on day 4",
pch="OD") +
scale_color_manual(values = c("red","blue","darkgreen")) +
theme(legend.position = "bottom",
legend.text = element_text(size = 30),
legend.title = element_text(size = 30),
axis.title = element_text(size = 30),
strip.background = element_blank(),
strip.text = element_text(size = 30, face = "italic"),
axis.text = element_text(size = 30))
Figure 1B - Interval plot of the proportion of flies
alive on day 4 post infection for various host genotypes (homozygous
serine dpt (wildtype, dptS69), homozygous arginine dpt (dptS69R), and
dpt null (Δdpt)), and variable pathogen strains (P.
burhodograneria Strain D, P. rettgeri, P.
sneebia, P. alcalifaciens) with different loads of
infective inoculum (OD).

Figure 2
Data_1$Treatment <- factor(Data_1$Treatment, levels = c("LB",
"P. alcalifaciens",
"P. burhodogranareia strD",
"P. rettgeri",
"P. sneebia"))
fit <- survfit(Surv(Day_Dead, censor) ~ OD_group + Treatment + Genotype,
data = Data_1)
# Figure 2A
Survival <- ggplot(data = surv_summary(fit),
aes(x = time,
y = surv,
group = interaction(Genotype,
OD_group))) +
facet_grid(~Treatment, scales = "free_x") +
coord_cartesian(ylim=c(0,1)) +
scale_x_discrete(
labels = function(x) {
if ("LB" %in% unique(Data$Treatment)) {
# For LB, show only '0'
ifelse(x == "0", "0", x)
} else {
# For other treatments, show the specific values
ifelse(x %in% c("0.01", "0.1", "0.5", "1", "4", "8"), x, NA)
}
}
) +
theme_classic() +
labs(x = "Days post infection",
y = "Survival") +
geom_step(aes(color = Genotype),
size = 1.5,
alpha=0.4) +
geom_point(aes(color = Genotype,
shape = OD_group),
size = 4) +
scale_shape_manual(values = c(4,1,2,0,5,7,8)) +
coord_cartesian(xlim = c(0, 7)) +
scale_color_manual(values = c("red",
"blue",
"darkgreen")) +
labs(pch = "OD") +
theme(legend.position = "bottom",
legend.text = element_text(size = 30),
legend.title = element_text(size = 30),
axis.title = element_text(size = 30),
strip.background = element_blank(),
strip.text = element_text(size = 30, face = "italic"),
axis.text = element_text(size = 30))
Fs$Treatment <- factor(Fs$Treatment, levels = c("LB",
"P. alcalifaciens",
"P. burhodogranareia strD",
"P. rettgeri",
"P. sneebia"))
# Figure 2B
FracSick <- ggplot(subset(Fs,
Day==4),
aes(x = as.factor(OD_group),
y = Fs,
pch = as.factor(OD_group),
col = Genotype)) +
stat_summary(geom = "errorbar",
position = position_dodge(0.7),
width = 0.3,
#aes(alpha = as.numeric(OD_group))
) +
stat_summary(geom = "point",
position = position_dodge(0.7),
size = 6,
#aes(alpha = as.numeric(OD_group))
) +
facet_grid(~Treatment, scales = "free_x") +
coord_cartesian(ylim=c(0,1)) +
scale_x_discrete(
labels = function(x) {
if ("LB" %in% unique(Data$Treatment)) {
# For LB, show only '0'
ifelse(x == "0", "0", x)
} else {
# For other treatments, show the specific values
ifelse(x %in% c("0.01", "0.1", "0.5", "1", "4", "8"), x, NA)
}
}
) +
scale_shape_manual(values = c(4,1,2,0,5,7,8)) +
theme_classic() +
labs(x = "OD",
y = "Fs on day 4",
pch="OD") +
scale_color_manual(values = c("red",
"blue",
"darkgreen")) +
theme(legend.position = "top",
legend.text = element_text(size = 30),
legend.title = element_text(size = 30),
axis.title = element_text(size = 30),
strip.background = element_blank(),
strip.text = element_text(size = 30, face = "italic"),
axis.text = element_text(size = 30)) +
guides(color=FALSE)
Data$Treatment <- factor(Data$Treatment, levels = c("LB",
"P. alcalifaciens",
"P. burhodogranareia strD",
"P. rettgeri",
"P. sneebia"))
# Figure 2C
Risk <- ggplot(Data,
aes(x = as.factor(OD_group),
y = Risk_Score,
pch = as.factor(OD_group),
col = Genotype)) +
stat_summary(geom = "errorbar",
position = position_dodge(0.7),
width = 0.3,
#aes(alpha = as.numeric(OD_group))
) +
stat_summary(geom = "point",
position = position_dodge(0.7),
size = 6,
#aes(alpha = as.numeric(OD_group))
) +
facet_grid(~Treatment, scales = "free_x") +
scale_x_discrete(
labels = function(x) {
if ("LB" %in% unique(Data$Treatment)) {
# For LB, show only '0'
ifelse(x == "0", "0", x)
} else {
# For other treatments, show the specific values
ifelse(x %in% c("0.01", "0.1", "0.5", "1", "4", "8"), x, NA)
}
}
) +
scale_shape_manual(values = c(4,1,2,0,5,7,8)) +
theme_classic() +
labs(x = "OD",
y = "Fs on day 4",
pch="OD") +
scale_color_manual(values = c("red",
"blue",
"darkgreen")) +
labs(x = "OD",
y = "Risk Score",
pch="OD") +
theme(legend.position = "top",
legend.text = element_text(size = 30),
legend.title = element_text(size = 30),
axis.title = element_text(size = 30),
strip.background = element_blank(),
strip.text = element_text(size = 30, face = "italic"),
axis.text = element_text(size = 30)) +
guides(color=FALSE)
Figure2 <- ggarrange(Survival, FracSick, Risk,
nrow = 3, ncol = 1,
common.legend = TRUE,
labels = "AUTO",
font.label = list(size = 40))
Figure 2 - Host-based measures of pathogenicity:
A. Survival curves,
B. Fraction symptomatic on Day4, and
C. Risk scores for infected flies of various host genotypes
(homozygous serine dpt (wildtype, dptS69), homozygous arginine dpt
(dptS69R), and dpt null (Δdpt)), and variable pathogen strains (P.
burhodograneria Strain D, P. rettgeri, P.
sneebia, P. alcalifaciens) with different loads of
infective inoculum (OD).

Figure 3
# Figure 3A
PP_D4 <- ggplot(Data,
aes(x = as.factor(OD_group),
y = PP_Day4,
pch = as.factor(OD_group),
col = Genotype)) +
stat_summary(geom = "errorbar",
position = position_dodge(0.7),
width = 0.3,
#aes(alpha = as.numeric(OD_group))
) +
stat_summary(geom = "point",
position = position_dodge(0.7),
size = 6,
#aes(alpha = as.numeric(OD_group))
) +
facet_grid(~Treatment, scales = "free_x") +
scale_x_discrete(
labels = function(x) {
if ("LB" %in% unique(Data$Treatment)) {
# For LB, show only '0'
ifelse(x == "0", "0", x)
} else {
# For other treatments, show the specific values
ifelse(x %in% c("0.01", "0.1", "0.5", "1", "4", "8"), x, NA)
}
}
) +
scale_shape_manual(values = c(4,1,2,0,5,7,8)) +
theme_classic() +
labs(x = "OD",
y = "Fs on day 4",
pch="OD") +
scale_color_manual(values = c("red",
"blue",
"darkgreen")) +
labs(x = "OD",
y = "PP on day 4",
pch="OD") +
theme(legend.position = "top",
legend.text = element_text(size = 30),
legend.title = element_text(size = 30),
axis.title = element_text(size = 30),
strip.background = element_blank(),
strip.text = element_text(size = 30, face = "italic"),
axis.text = element_text(size = 30))
# Figure 3B
PP_T <- ggplot(Data,
aes(x = as.factor(OD_group),
y = PP_T,
pch = as.factor(OD_group),
col = Genotype)) +
stat_summary(geom = "errorbar",
position = position_dodge(0.7),
width = 0.3,
#aes(alpha = as.numeric(OD_group))
) +
stat_summary(geom = "point",
position = position_dodge(0.7),
size = 6,
#aes(alpha = as.numeric(OD_group))
) +
facet_grid(~Treatment, scales = "free_x") +
scale_x_discrete(
labels = function(x) {
if ("LB" %in% unique(Data$Treatment)) {
# For LB, show only '0'
ifelse(x == "0", "0", x)
} else {
# For other treatments, show the specific values
ifelse(x %in% c("0.01", "0.1", "0.5", "1", "4", "8"), x, NA)
}
}
) +
scale_shape_manual(values = c(4,1,2,0,5,7,8)) +
theme_classic() +
labs(x = "OD",
y = "Fs on day 4",
pch="OD") +
scale_color_manual(values = c("red",
"blue",
"darkgreen")) +
labs(x = "OD",
y = "PP_T",
pch="OD") +
theme(legend.position = "top",
legend.text = element_text(size = 30),
legend.title = element_text(size = 30),
axis.title = element_text(size = 30),
strip.background = element_blank(),
strip.text = element_text(size = 30, face = "italic"),
axis.text = element_text(size = 30))
# Figure 3C
Fs_T <- ggplot(Data,
aes(x = as.factor(OD_group),
y = Fs_T,
pch = as.factor(OD_group),
col = Genotype)) +
stat_summary(geom = "errorbar",
position = position_dodge(0.7),
width = 0.3,
#aes(alpha = as.numeric(OD_group))
) +
stat_summary(geom = "point",
position = position_dodge(0.7),
size = 6,
#aes(alpha = as.numeric(OD_group))
) +
facet_grid(~Treatment, scales = "free_x") +
coord_cartesian(ylim=c(0,1)) +
scale_x_discrete(
labels = function(x) {
if ("LB" %in% unique(Data$Treatment)) {
# For LB, show only '0'
ifelse(x == "0", "0", x)
} else {
# For other treatments, show the specific values
ifelse(x %in% c("0.01", "0.1", "0.5", "1", "4", "8"), x, NA)
}
}
) +
scale_shape_manual(values = c(4,1,2,0,5,7,8)) +
theme_classic() +
labs(x = "OD",
y = "Fs on day 4",
pch="OD") +
scale_color_manual(values = c("red",
"blue",
"darkgreen")) +
labs(x = "OD",
y = "Fs_T",
pch="OD") +
theme(legend.position = "top",
legend.text = element_text(size = 30),
legend.title = element_text(size = 30),
axis.title = element_text(size = 30),
strip.background = element_blank(),
strip.text = element_text(size = 30, face = "italic"),
axis.text = element_text(size = 30))
Figure3 <- ggarrange(PP_D4, PP_T, Fs_T,
nrow = 3, ncol = 1,
common.legend = TRUE,
labels = "AUTO",
font.label = list(size = 40))
Figure 3
- A. Pathogen potential on Day 4 post-infection,
- B. Pathogen potential on median survival day (T), and
- C. Fraction symptomatic on median survival day (T) for infected
flies of various host genotypes (homozygous serine dpt (wildtype,
dptS69), homozygous arginine dpt (dptS69R), and dpt null (Δdpt)), and
variable pathogen strains (P. burhodograneria Strain D, P.
rettgeri, P. sneebia, P. alcalifaciens) with
different loads of infective inoculum (OD).

Figure 4
Figure 4A
# Variables of interest
variables_of_interest <- c("Risk_Score",
"Prop_Day4_Alive",
"Fs_Day4",
"PP_Day4",
"Fs_T",
"PP_T")
# Create an empty dataframe to store p-values matrices
Data_stats_1 <- data.frame(
variable = character(),
OD_group = character(),
Genotype = character(),
Treatment = character(),
`OD_group*Genotype` = character(),
`OD_group*Treatment` = character(),
`Genotype*Treatment` = character(),
`OD_group*Genotype*Treatment` = character(),
stringsAsFactors = FALSE
)
TableS3 <- list()
# Loop through each variable of interest
for (variable in variables_of_interest) {
# Fit an ANOVA model
model <- aov(as.formula(paste(variable, "~ OD_group * Genotype * Treatment + Error(Vial_ID+Date)")), data = Data)
# Extract ANOVA table
anova_table <- as.data.frame(summary(model$Within)[[1]])
TableS3[[variable]] <- anova_table
# Create a new row to append
new_row <- data.frame(
variable = variable,
OD_group = anova_table$`Pr(>F)`[1],
Genotype = anova_table$`Pr(>F)`[2],
Treatment = anova_table$`Pr(>F)`[3],
`OD_group*Genotype` = anova_table$`Pr(>F)`[4],
`OD_group*Treatment` = anova_table$`Pr(>F)`[5],
`Genotype*Treatment` = anova_table$`Pr(>F)`[6],
`OD_group*Genotype*Treatment` = anova_table$`Pr(>F)`[7]
)
# Append the new row to Data_stats_1
Data_stats_1 <- rbind(Data_stats_1, new_row)
}
Risk_Score_anova <- TableS3[["Risk_Score"]]
Prop_Day4_Alive_anova <- TableS3[["Prop_Day4_Alive"]]
Fs_Day4_anova <- TableS3[["Fs_Day4"]]
PP_Day4_anova <- TableS3[["PP_Day4"]]
Fs_T_anova <- TableS3[["Fs_T"]]
PP_T_anova <- TableS3[["PP_T"]]
Data_stats <- as.data.frame(melt(Data_stats_1, id.vars = c("variable"), variable.name = "vars"))
Data_stats$pval <- ifelse(Data_stats$value>0.05,"Not Significant",
ifelse(Data_stats$value<0.05 &
Data_stats$value>0.01, "P<0.05",
ifelse(Data_stats$value<0.01 &
Data_stats$value>0.001, "P<0.01",
ifelse(Data_stats$value<0.001, "P<0.001","NA"))))
custom_labels_vars <- c("OD_group" = "OD",
"Genotype" = "Host genotype",
"Treatment" = "Pathogen species",
"OD_group.Genotype" = "OD * Host genotype",
"OD_group.Treatment" = "OD * Pathogen species",
"Genotype.Treatment" = "Host genotype * Pathogen species",
"OD_group.Genotype.Treatment" = "OD * Host genotype * Pathogen species")
custom_labels_variable <- c("Fs_Day4" = "Fs on Day4",
"Fs_T" = "Fs_T",
"PP_Day4" = "PP on Day4",
"PP_T" = "PP_T",
"Prop_Day4_Alive" = "Proportion alive on Day4",
"Risk_Score" = "Risk score")
Fig4A <- ggplot(Data_stats, aes(x=as.factor(vars),
y=as.factor(variable),
fill=as.factor(pval))) +
geom_tile() +
theme_bw() +
labs(x="Factor",
y="Measure of pathogenicity",
fill="p-value") +
scale_x_discrete(labels = custom_labels_vars) +
scale_y_discrete(labels = custom_labels_variable) +
scale_fill_manual(values = c("white","maroon","red","pink")) +
theme(axis.text.y = element_text(size=40,
colour = c("darkorchid2","darkorchid2","orange","orange","darkorchid2","darkorchid2")),
axis.text.x = element_text(size=40, angle = 45, vjust = 1, hjust = 1,
colour = c("black","darkorchid2","orange","black","black","black","black")),
axis.title = element_text(size = 40),
legend.position = "top",
legend.text = element_text(size = 40),
legend.title = element_text(size = 40))
Figure 4B
# Select the columns for analysis
Corrdata <- Data[,c(27,34,42,46,48,49)]
# Compute the correlation matrix
M <- cor(Corrdata)
TableS4 <- as.data.frame(M)
# Compute the R^2 matrix
R2 <- M^2
# Melt R2 matrix into a long format
df <- melt(R2)
# Filter out the upper triangle of the matrix
df <- df[as.numeric(df$Var1) >= as.numeric(df$Var2), ]
# Correlation plot - Figure 4B
Fig4B <- ggplot(data = df,
aes(Var1, Var2,
fill = value)) +
geom_tile(color = "black") +
scale_fill_gradient2(low = "white",
high = "lightblue",
mid = "grey90",
midpoint = 0.5) +
geom_text(aes(label = round(value, 2)),
size = 20) +
theme_classic() +
theme(axis.text.x = element_text(colour = c("darkorchid2",
"orange",
"darkorchid2",
"darkorchid2",
"darkorchid2",
"orange"),
size = 40,
angle = 45,
vjust = 1,
hjust = 1,),
axis.text.y = element_text(colour = c("darkorchid2",
"orange",
"darkorchid2",
"darkorchid2",
"darkorchid2",
"orange"),
size = 40),
legend.position = "none") +
scale_x_discrete(labels = custom_labels_variable) +
scale_y_discrete(labels = custom_labels_variable) +
labs(x="",
y="")
Figure4 <- ggarrange(Fig4A, Fig4B,
nrow = 1, ncol = 2,
labels = "AUTO",
font.label = list(size = 50))
Figure 4A P-values from ANOVA to indicate the
significance of each factor corresponding to each measure of
pathogenicity
Figure 4B Heatmap showing R2 (R=correlation
coefficient) values to assess the differences between different measures
of host susceptibility (survival, risk scores, and fraction symptomatic)
and pathogen virulence (pathogen potential)
[Prop_Day4_Alive: Proportion survival on Day 4 post-infection,
PP_Day4: Pathogen potential on Day 4, PP_T: Pathogen potential on median
survival day, Fs_Day4: Fraction symptomatic on Day 4, Fs_T: Fraction
symptomatic on median survival day].

Statistical Analysis
Table S1
Table S1: Statistical analysis for analyzing the
effect of OD and treatment (pathogen species) on log-transformed
bacterial load. Data was fitted to a linear mixed-effects model.
[Model: lm( log(Bacterial_Load_CFU/fly+1) ~ OD +
Treatment)]
Model1 <- lm(`Bacterial_Load_CFU/fly`~OD+Treatment,data = LoadData)
MSum1 <- summary(Model1)
TableS1 <- as.data.frame(MSum1$coefficients)
print(TableS1)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -173.1791 1319.6811 -0.1312280 0.896057039
## OD 976.3464 281.7327 3.4655067 0.001012516
## TreatmentP. sneebia 191.8333 1632.8846 0.1174813 0.906891692
Table S2
Table S2: Statistical analysis for analyzing the
effect of OD, treatment (pathogen species), and host genotype on
survival. Data was fitted to a mixed-effects cox-proportional hazards
model.
[Model: coxph(Survival ~ OD * Treatment * Genotype]
cox <- coxph(Surv(Day_Dead, censor) ~ OD_group * Treatment * Genotype +
frailty(Date),
data = Data_1)
cox_summary <- summary(cox)
TableS2 <- as.data.frame(cox_summary$coefficients)
print(TableS2)
## coef se(coef) se2
## OD_group0.01 2.465579129 4.924019e-01 4.882623e-01
## OD_group0.1 2.543120127 4.850304e-01 4.808286e-01
## OD_group0.5 2.715388454 4.744765e-01 4.744764e-01
## OD_group1 -0.141432999 4.765642e-01 4.765642e-01
## OD_group4 2.853860034 4.698408e-01 4.698407e-01
## OD_group8 4.149345144 4.736175e-01 4.693276e-01
## TreatmentP..alcalifaciens -1.671417849 6.982090e+05 6.982090e+05
## TreatmentP..burhodogranar -1.914697543 2.427268e-01 2.328724e-01
## TreatmentP..rettgeri -1.246073645 3.722759e+05 3.722759e+05
## TreatmentP..sneebia NA 0.000000e+00 0.000000e+00
## GenotypeDptR -0.382912516 6.732125e-01 6.712090e-01
## GenotypeDptS -0.309707667 6.717803e-01 6.709762e-01
## frailty.Date. NA NA NA
## OD_group0.01.TreatmentP.. 2.399820178 6.982090e+05 6.982090e+05
## OD_group0.1.TreatmentP..a 2.156073003 6.982090e+05 6.982090e+05
## OD_group0.5.TreatmentP..a 2.165868733 6.982090e+05 6.982090e+05
## OD_group1.TreatmentP..alc 5.073849259 6.982090e+05 6.982090e+05
## OD_group4.TreatmentP..alc 2.598092750 6.982090e+05 6.982090e+05
## OD_group8.TreatmentP..alc 1.350307927 6.982090e+05 6.982090e+05
## OD_group0.01.TreatmentP...1 1.371817724 3.648377e-01 3.648376e-01
## OD_group0.1.TreatmentP..b 0.752297228 3.760251e-01 3.760251e-01
## OD_group0.5.TreatmentP..b 0.541166546 3.772288e-01 3.718072e-01
## OD_group1.TreatmentP..bur 4.213991301 3.482430e-01 3.423709e-01
## OD_group4.TreatmentP..bur 0.996788604 3.425627e-01 3.365727e-01
## OD_group8.TreatmentP..bur NA 0.000000e+00 0.000000e+00
## OD_group0.01.TreatmentP...2 0.562711101 3.722759e+05 3.722759e+05
## OD_group0.1.TreatmentP..r 1.452804529 3.722759e+05 3.722759e+05
## OD_group0.5.TreatmentP..r 1.502594030 3.722759e+05 3.722759e+05
## OD_group1.TreatmentP..ret 4.785627270 3.722759e+05 3.722759e+05
## OD_group4.TreatmentP..ret 1.264499761 3.722759e+05 3.722759e+05
## OD_group8.TreatmentP..ret 0.171801966 3.722759e+05 3.722759e+05
## OD_group0.01.TreatmentP...3 NA 0.000000e+00 0.000000e+00
## OD_group0.1.TreatmentP..s NA 0.000000e+00 0.000000e+00
## OD_group0.5.TreatmentP..s NA 0.000000e+00 0.000000e+00
## OD_group1.TreatmentP..sne NA 0.000000e+00 0.000000e+00
## OD_group4.TreatmentP..sne NA 0.000000e+00 0.000000e+00
## OD_group8.TreatmentP..sne NA 0.000000e+00 0.000000e+00
## OD_group0.01.GenotypeDptR 0.470713354 7.281076e-01 7.262556e-01
## OD_group0.1.GenotypeDptR 0.883859006 7.142435e-01 7.123480e-01
## OD_group0.5.GenotypeDptR 1.055857403 7.045312e-01 7.045311e-01
## OD_group1.GenotypeDptR 0.937833652 7.075987e-01 7.070637e-01
## OD_group4.GenotypeDptR 0.452329473 7.009321e-01 7.003985e-01
## OD_group8.GenotypeDptR 0.420449012 7.025855e-01 7.006654e-01
## OD_group0.01.GenotypeDptS 0.449644062 7.240409e-01 7.232943e-01
## OD_group0.1.GenotypeDptS 0.266941237 7.143786e-01 7.136220e-01
## OD_group0.5.GenotypeDptS 0.882105746 7.048471e-01 7.048330e-01
## OD_group1.GenotypeDptS 1.035919516 7.059839e-01 7.059634e-01
## OD_group4.GenotypeDptS 0.121836449 7.015740e-01 7.015551e-01
## OD_group8.GenotypeDptS -0.060228758 7.018380e-01 7.010733e-01
## TreatmentP..alcalifaciens.1 -0.085236782 2.835734e-01 2.835734e-01
## TreatmentP..burhodogranar.1 -0.042165244 3.289776e-01 3.289776e-01
## TreatmentP..rettgeri.Geno 0.004937081 2.843097e-01 2.843097e-01
## TreatmentP..sneebia.Genot NA 0.000000e+00 0.000000e+00
## TreatmentP..alcalifaciens.2 0.322236139 2.851122e-01 2.851122e-01
## TreatmentP..burhodogranar.2 0.276863162 3.338959e-01 3.338959e-01
## TreatmentP..rettgeri.Geno.1 -0.187161768 2.937825e-01 2.937825e-01
## TreatmentP..sneebia.Genot.1 NA 0.000000e+00 0.000000e+00
## OD_group0.01.TreatmentP...4 -0.314937199 4.504374e-01 4.504374e-01
## OD_group0.1.TreatmentP..a.1 -0.330612917 4.257572e-01 4.257572e-01
## OD_group0.5.TreatmentP..a.1 -0.106091002 4.127077e-01 4.094277e-01
## OD_group1.TreatmentP..alc.1 0.097547468 4.145128e-01 4.137762e-01
## OD_group4.TreatmentP..alc.1 -0.029708361 4.019184e-01 4.011670e-01
## OD_group8.TreatmentP..alc.1 NA 0.000000e+00 0.000000e+00
## OD_group0.01.TreatmentP...5 -0.710053077 5.363291e-01 5.363291e-01
## OD_group0.1.TreatmentP..b.1 0.158952671 5.076546e-01 5.076546e-01
## OD_group0.5.TreatmentP..b.1 -0.847586508 5.366671e-01 5.341718e-01
## OD_group1.TreatmentP..bur.1 -2.066817617 5.528053e-01 5.522705e-01
## OD_group4.TreatmentP..bur.1 0.259844109 4.714907e-01 4.708480e-01
## OD_group8.TreatmentP..bur.1 NA 0.000000e+00 0.000000e+00
## OD_group0.01.TreatmentP...6 0.224937309 4.867779e-01 4.867779e-01
## OD_group0.1.TreatmentP..r.1 -0.304239563 4.286863e-01 4.286863e-01
## OD_group0.5.TreatmentP..r.1 -0.565978533 4.137873e-01 4.105542e-01
## OD_group1.TreatmentP..ret.1 -0.677703633 4.150179e-01 4.143087e-01
## OD_group4.TreatmentP..ret.1 -0.272944117 4.111758e-01 4.104495e-01
## OD_group8.TreatmentP..ret.1 NA 0.000000e+00 0.000000e+00
## OD_group0.01.TreatmentP...7 NA 0.000000e+00 0.000000e+00
## OD_group0.1.TreatmentP..s.1 NA 0.000000e+00 0.000000e+00
## OD_group0.5.TreatmentP..s.1 NA 0.000000e+00 0.000000e+00
## OD_group1.TreatmentP..sne.1 NA 0.000000e+00 0.000000e+00
## OD_group4.TreatmentP..sne.1 NA 0.000000e+00 0.000000e+00
## OD_group8.TreatmentP..sne.1 NA 0.000000e+00 0.000000e+00
## OD_group0.01.TreatmentP...8 -1.431941830 4.582501e-01 4.525892e-01
## OD_group0.1.TreatmentP..a.2 -0.537402197 4.368553e-01 4.309135e-01
## OD_group0.5.TreatmentP..a.2 -0.826002197 4.170152e-01 4.098035e-01
## OD_group1.TreatmentP..alc.2 -1.058159287 4.204203e-01 4.133327e-01
## OD_group4.TreatmentP..alc.2 -0.265166412 4.054168e-01 4.044784e-01
## OD_group8.TreatmentP..alc.2 NA 0.000000e+00 0.000000e+00
## OD_group0.01.TreatmentP...9 -1.533479745 5.548007e-01 5.548007e-01
## OD_group0.1.TreatmentP..b.2 -0.542555529 5.454616e-01 5.454616e-01
## OD_group0.5.TreatmentP..b.2 -1.013040316 5.336373e-01 5.328856e-01
## OD_group1.TreatmentP..bur.2 -2.096057726 5.177448e-01 5.170239e-01
## OD_group4.TreatmentP..bur.2 0.043710955 4.796295e-01 4.788321e-01
## OD_group8.TreatmentP..bur.2 NA 0.000000e+00 0.000000e+00
## OD_group0.01.TreatmentP...10 0.030270044 4.943688e-01 4.943688e-01
## OD_group0.1.TreatmentP..r.2 -0.111723849 4.428810e-01 4.428810e-01
## OD_group0.5.TreatmentP..r.2 -0.646098913 4.222030e-01 4.212477e-01
## OD_group1.TreatmentP..ret.2 -0.975086453 4.206219e-01 4.197294e-01
## OD_group4.TreatmentP..ret.2 -0.366886886 4.288274e-01 4.279366e-01
## OD_group8.TreatmentP..ret.2 NA 0.000000e+00 0.000000e+00
## OD_group0.01.TreatmentP...11 NA 0.000000e+00 0.000000e+00
## OD_group0.1.TreatmentP..s.2 NA 0.000000e+00 0.000000e+00
## OD_group0.5.TreatmentP..s.2 NA 0.000000e+00 0.000000e+00
## OD_group1.TreatmentP..sne.2 NA 0.000000e+00 0.000000e+00
## OD_group4.TreatmentP..sne.2 NA 0.000000e+00 0.000000e+00
## OD_group8.TreatmentP..sne.2 NA 0.000000e+00 0.000000e+00
## Chisq DF p
## OD_group0.01 2.507255e+01 1.00000000 5.521322e-07
## OD_group0.1 2.749134e+01 1.00000000 1.577995e-07
## OD_group0.5 3.275175e+01 1.00000000 1.047119e-08
## OD_group1 8.807622e-02 1.00000000 7.666372e-01
## OD_group4 3.689469e+01 1.00000000 1.246852e-09
## OD_group8 7.675448e+01 1.00000000 1.935887e-18
## TreatmentP..alcalifaciens 5.730588e-12 1.00000000 9.999981e-01
## TreatmentP..burhodogranar 6.222500e+01 1.00000000 3.063741e-15
## TreatmentP..rettgeri 1.120360e-11 1.00000000 9.999973e-01
## TreatmentP..sneebia NA 1.00000000 NA
## GenotypeDptR 3.235153e-01 1.00000000 5.695027e-01
## GenotypeDptS 2.125442e-01 1.00000000 6.447805e-01
## frailty.Date. 3.525123e+00 -0.08584406 2.353071e-02
## OD_group0.01.TreatmentP.. 1.181372e-11 1.00000000 9.999973e-01
## OD_group0.1.TreatmentP..a 9.535776e-12 1.00000000 9.999975e-01
## OD_group0.5.TreatmentP..a 9.622621e-12 1.00000000 9.999975e-01
## OD_group1.TreatmentP..alc 5.280855e-11 1.00000000 9.999942e-01
## OD_group4.TreatmentP..alc 1.384645e-11 1.00000000 9.999970e-01
## OD_group8.TreatmentP..alc 3.740199e-12 1.00000000 9.999985e-01
## OD_group0.01.TreatmentP...1 1.413818e+01 1.00000000 1.698613e-04
## OD_group0.1.TreatmentP..b 4.002628e+00 1.00000000 4.542937e-02
## OD_group0.5.TreatmentP..b 2.058032e+00 1.00000000 1.514056e-01
## OD_group1.TreatmentP..bur 1.464275e+02 1.00000000 1.046858e-33
## OD_group4.TreatmentP..bur 8.466930e+00 1.00000000 3.616612e-03
## OD_group8.TreatmentP..bur NA 1.00000000 NA
## OD_group0.01.TreatmentP...2 2.284763e-12 1.00000000 9.999988e-01
## OD_group0.1.TreatmentP..r 1.522946e-11 1.00000000 9.999969e-01
## OD_group0.5.TreatmentP..r 1.629121e-11 1.00000000 9.999968e-01
## OD_group1.TreatmentP..ret 1.652524e-10 1.00000000 9.999897e-01
## OD_group4.TreatmentP..ret 1.153739e-11 1.00000000 9.999973e-01
## OD_group8.TreatmentP..ret 2.129739e-13 1.00000000 9.999996e-01
## OD_group0.01.TreatmentP...3 NA 1.00000000 NA
## OD_group0.1.TreatmentP..s NA 1.00000000 NA
## OD_group0.5.TreatmentP..s NA 1.00000000 NA
## OD_group1.TreatmentP..sne NA 1.00000000 NA
## OD_group4.TreatmentP..sne NA 1.00000000 NA
## OD_group8.TreatmentP..sne NA 1.00000000 NA
## OD_group0.01.GenotypeDptR 4.179477e-01 1.00000000 5.179629e-01
## OD_group0.1.GenotypeDptR 1.531346e+00 1.00000000 2.159105e-01
## OD_group0.5.GenotypeDptR 2.246002e+00 1.00000000 1.339601e-01
## OD_group1.GenotypeDptR 1.756619e+00 1.00000000 1.850468e-01
## OD_group4.GenotypeDptR 4.164453e-01 1.00000000 5.187161e-01
## OD_group8.GenotypeDptR 3.581198e-01 1.00000000 5.495523e-01
## OD_group0.01.GenotypeDptS 3.856661e-01 1.00000000 5.345863e-01
## OD_group0.1.GenotypeDptS 1.396286e-01 1.00000000 7.086505e-01
## OD_group0.5.GenotypeDptS 1.566215e+00 1.00000000 2.107575e-01
## OD_group1.GenotypeDptS 2.153091e+00 1.00000000 1.422831e-01
## OD_group4.GenotypeDptS 3.015835e-02 1.00000000 8.621314e-01
## OD_group8.GenotypeDptS 7.364344e-03 1.00000000 9.316129e-01
## TreatmentP..alcalifaciens.1 9.034895e-02 1.00000000 7.637340e-01
## TreatmentP..burhodogranar.1 1.642770e-02 1.00000000 8.980140e-01
## TreatmentP..rettgeri.Geno 3.015486e-04 1.00000000 9.861453e-01
## TreatmentP..sneebia.Genot NA 1.00000000 NA
## TreatmentP..alcalifaciens.2 1.277370e+00 1.00000000 2.583886e-01
## TreatmentP..burhodogranar.2 6.875560e-01 1.00000000 4.069969e-01
## TreatmentP..rettgeri.Geno.1 4.058658e-01 1.00000000 5.240754e-01
## TreatmentP..sneebia.Genot.1 NA 1.00000000 NA
## OD_group0.01.TreatmentP...4 4.888538e-01 1.00000000 4.844391e-01
## OD_group0.1.TreatmentP..a.1 6.029975e-01 1.00000000 4.374366e-01
## OD_group0.5.TreatmentP..a.1 6.608028e-02 1.00000000 7.971318e-01
## OD_group1.TreatmentP..alc.1 5.538041e-02 1.00000000 8.139523e-01
## OD_group4.TreatmentP..alc.1 5.463633e-03 1.00000000 9.410769e-01
## OD_group8.TreatmentP..alc.1 NA 1.00000000 NA
## OD_group0.01.TreatmentP...5 1.752746e+00 1.00000000 1.855319e-01
## OD_group0.1.TreatmentP..b.1 9.803904e-02 1.00000000 7.541957e-01
## OD_group0.5.TreatmentP..b.1 2.494354e+00 1.00000000 1.142552e-01
## OD_group1.TreatmentP..bur.1 1.397848e+01 1.00000000 1.849152e-04
## OD_group4.TreatmentP..bur.1 3.037243e-01 1.00000000 5.815570e-01
## OD_group8.TreatmentP..bur.1 NA 1.00000000 NA
## OD_group0.01.TreatmentP...6 2.135312e-01 1.00000000 6.440137e-01
## OD_group0.1.TreatmentP..r.1 5.036770e-01 1.00000000 4.778889e-01
## OD_group0.5.TreatmentP..r.1 1.870878e+00 1.00000000 1.713745e-01
## OD_group1.TreatmentP..ret.1 2.666527e+00 1.00000000 1.024794e-01
## OD_group4.TreatmentP..ret.1 4.406487e-01 1.00000000 5.068095e-01
## OD_group8.TreatmentP..ret.1 NA 1.00000000 NA
## OD_group0.01.TreatmentP...7 NA 1.00000000 NA
## OD_group0.1.TreatmentP..s.1 NA 1.00000000 NA
## OD_group0.5.TreatmentP..s.1 NA 1.00000000 NA
## OD_group1.TreatmentP..sne.1 NA 1.00000000 NA
## OD_group4.TreatmentP..sne.1 NA 1.00000000 NA
## OD_group8.TreatmentP..sne.1 NA 1.00000000 NA
## OD_group0.01.TreatmentP...8 9.764399e+00 1.00000000 1.779236e-03
## OD_group0.1.TreatmentP..a.2 1.513295e+00 1.00000000 2.186370e-01
## OD_group0.5.TreatmentP..a.2 3.923364e+00 1.00000000 4.761949e-02
## OD_group1.TreatmentP..alc.2 6.334828e+00 1.00000000 1.183896e-02
## OD_group4.TreatmentP..alc.2 4.277930e-01 1.00000000 5.130739e-01
## OD_group8.TreatmentP..alc.2 NA 1.00000000 NA
## OD_group0.01.TreatmentP...9 7.639801e+00 1.00000000 5.709422e-03
## OD_group0.1.TreatmentP..b.2 9.893728e-01 1.00000000 3.198957e-01
## OD_group0.5.TreatmentP..b.2 3.603803e+00 1.00000000 5.764755e-02
## OD_group1.TreatmentP..bur.2 1.638985e+01 1.00000000 5.156051e-05
## OD_group4.TreatmentP..bur.2 8.305557e-03 1.00000000 9.273855e-01
## OD_group8.TreatmentP..bur.2 NA 1.00000000 NA
## OD_group0.01.TreatmentP...10 3.749074e-03 1.00000000 9.511763e-01
## OD_group0.1.TreatmentP..r.2 6.363818e-02 1.00000000 8.008354e-01
## OD_group0.5.TreatmentP..r.2 2.341831e+00 1.00000000 1.259415e-01
## OD_group1.TreatmentP..ret.2 5.374061e+00 1.00000000 2.043834e-02
## OD_group4.TreatmentP..ret.2 7.319801e-01 1.00000000 3.922421e-01
## OD_group8.TreatmentP..ret.2 NA 1.00000000 NA
## OD_group0.01.TreatmentP...11 NA 1.00000000 NA
## OD_group0.1.TreatmentP..s.2 NA 1.00000000 NA
## OD_group0.5.TreatmentP..s.2 NA 1.00000000 NA
## OD_group1.TreatmentP..sne.2 NA 1.00000000 NA
## OD_group4.TreatmentP..sne.2 NA 1.00000000 NA
## OD_group8.TreatmentP..sne.2 NA 1.00000000 NA
Table S3
Table S3: ANOVA for analyzing the effect of OD,
treatment (pathogen species), and host genotype on different measures of
pathogenicity, corresponding to Figure 4A.
[Model: aov(“Measure of pathogenicity” ~ OD * Treatment *
Genotype]
# Variables of interest
variables_of_interest <- c("Risk_Score",
"Prop_Day4_Alive",
"Fs_Day4",
"PP_Day4",
"Fs_T",
"PP_T")
# Create an empty dataframe to store p-values matrices
Data_stats_1 <- data.frame(
variable = character(),
OD_group = character(),
Genotype = character(),
Treatment = character(),
`OD_group*Genotype` = character(),
`OD_group*Treatment` = character(),
`Genotype*Treatment` = character(),
`OD_group*Genotype*Treatment` = character(),
stringsAsFactors = FALSE
)
TableS3 <- list()
# Loop through each variable of interest
for (variable in variables_of_interest) {
# Fit an ANOVA model
model <- aov(as.formula(paste(variable, "~ OD_group * Genotype * Treatment + Error(Vial_ID+Date)")), data = Data)
# Extract ANOVA table
anova_table <- as.data.frame(summary(model$Within)[[1]])
TableS3[[variable]] <- anova_table
# Create a new row to append
new_row <- data.frame(
variable = variable,
OD_group = anova_table$`Pr(>F)`[1],
Genotype = anova_table$`Pr(>F)`[2],
Treatment = anova_table$`Pr(>F)`[3],
`OD_group*Genotype` = anova_table$`Pr(>F)`[4],
`OD_group*Treatment` = anova_table$`Pr(>F)`[5],
`Genotype*Treatment` = anova_table$`Pr(>F)`[6],
`OD_group*Genotype*Treatment` = anova_table$`Pr(>F)`[7]
)
# Append the new row to Data_stats_1
Data_stats_1 <- rbind(Data_stats_1, new_row)
}
Risk_Score_anova <- TableS3[["Risk_Score"]]
Prop_Day4_Alive_anova <- TableS3[["Prop_Day4_Alive"]]
Fs_Day4_anova <- TableS3[["Fs_Day4"]]
PP_Day4_anova <- TableS3[["PP_Day4"]]
Fs_T_anova <- TableS3[["Fs_T"]]
PP_T_anova <- TableS3[["PP_T"]]
A. Risk score
## Df Sum Sq Mean Sq F value
## OD_group 1 1.721013e+17 1.721013e+17 188.7850978
## Genotype 2 2.380668e+16 1.190334e+16 13.0572697
## Treatment 1 8.631467e+16 8.631467e+16 94.6821660
## OD_group:Genotype 2 8.151767e+14 4.075883e+14 0.4471007
## OD_group:Treatment 3 3.397697e+16 1.132566e+16 12.4235869
## Genotype:Treatment 8 2.442590e+16 3.053237e+15 3.3492234
## OD_group:Genotype:Treatment 6 2.156488e+16 3.594146e+15 3.9425694
## Residuals 343 3.126875e+17 9.116254e+14 NA
## Pr(>F)
## OD_group 1.569612e-34
## Genotype 3.426131e-06
## Treatment 6.399903e-20
## OD_group:Genotype 6.398517e-01
## OD_group:Treatment 9.870310e-08
## Genotype:Treatment 1.041046e-03
## OD_group:Genotype:Treatment 7.938854e-04
## Residuals NA
B. Proportion alive on day 4
## Df Sum Sq Mean Sq F value Pr(>F)
## OD_group 1 2.84511946 2.84511946 146.4453151 2.580869e-28
## Genotype 2 0.22638178 0.11319089 5.8262143 3.248910e-03
## Treatment 1 6.24734970 6.24734970 321.5664960 3.371913e-51
## OD_group:Genotype 2 0.00645126 0.00322563 0.1660311 8.470879e-01
## OD_group:Treatment 3 0.36332099 0.12110700 6.2336758 3.933111e-04
## Genotype:Treatment 8 0.27880303 0.03485038 1.7938350 7.717203e-02
## OD_group:Genotype:Treatment 6 0.17242848 0.02873808 1.4792198 1.843522e-01
## Residuals 343 6.66375687 0.01942786 NA NA
D. Pathogen potential on day 4 (PP_Day4)
## Df Sum Sq Mean Sq F value
## OD_group 1 419.1892066 419.18920656 236.74644014
## Genotype 2 1.7895488 0.89477438 0.50534376
## Treatment 1 72.7881222 72.78812220 41.10871307
## OD_group:Genotype 2 0.1078957 0.05394785 0.03046825
## OD_group:Treatment 3 68.1946760 22.73155868 12.83815402
## Genotype:Treatment 8 27.7610991 3.47013739 1.95983737
## OD_group:Genotype:Treatment 6 28.9251722 4.82086203 2.72268920
## Residuals 343 607.3244344 1.77062517 NA
## Pr(>F)
## OD_group 5.423226e-41
## Genotype 6.037466e-01
## Treatment 4.768189e-10
## OD_group:Genotype 9.699939e-01
## OD_group:Treatment 5.717453e-08
## Genotype:Treatment 5.072378e-02
## OD_group:Genotype:Treatment 1.347991e-02
## Residuals NA
F. Fraction symptomatic on day 4 (Fs_Day4)
## Df Sum Sq Mean Sq F value Pr(>F)
## OD_group 1 0.39252080 0.392520798 29.3149653 1.158839e-07
## Genotype 2 0.05260709 0.026303544 1.9644500 1.418076e-01
## Treatment 1 0.98637215 0.986372151 73.6660720 3.264645e-16
## OD_group:Genotype 2 0.01769875 0.008849373 0.6609053 5.170400e-01
## OD_group:Treatment 3 0.05409545 0.018031818 1.3466856 2.590671e-01
## Genotype:Treatment 8 0.24379423 0.030474279 2.2759366 2.202909e-02
## OD_group:Genotype:Treatment 6 0.16760849 0.027934748 2.0862746 5.428488e-02
## Residuals 343 4.59269293 0.013389775 NA NA
Table S4
Table S4: Correlation matrix for different measures
of pathogenicity containing R values corresponding to Figure 4B.
# Select the columns for analysis
Corrdata <- Data[,c(27,34,42,46,48:50)]
# Compute the correlation matrix
M <- cor(Corrdata)
TableS4 <- as.data.frame(M)
print(TableS4)
## Fs_Day4 PP_Day4 Prop_Day4_Alive Risk_Score Fs_T
## Fs_Day4 1.0000000 0.3035784 -0.7881867 0.6467606 0.6756411
## PP_Day4 0.3035784 1.0000000 -0.3719103 0.2167741 0.2059444
## Prop_Day4_Alive -0.7881867 -0.3719103 1.0000000 -0.8139554 -0.5490079
## Risk_Score 0.6467606 0.2167741 -0.8139554 1.0000000 0.6024717
## Fs_T 0.6756411 0.2059444 -0.5490079 0.6024717 1.0000000
## PP_T 0.1941431 0.8643599 -0.2450816 0.1913188 0.3243466
## Fs_T/T 0.6789812 0.3029858 -0.8681956 0.8430373 0.5795366
## PP_T Fs_T/T
## Fs_Day4 0.1941431 0.6789812
## PP_Day4 0.8643599 0.3029858
## Prop_Day4_Alive -0.2450816 -0.8681956
## Risk_Score 0.1913188 0.8430373
## Fs_T 0.3243466 0.5795366
## PP_T 1.0000000 0.2066825
## Fs_T/T 0.2066825 1.0000000